Analysis of blood methylation quantitative trait loci in East Asians reveals ancestry-specific impacts on complex traits.

Methylation quantitative trait loci (mQTLs) are essential for understanding the role of DNA methylation changes in genetic predisposition, yet they have not been fully characterized in East Asians (EAs). Here we identified mQTLs in whole blood from 3,523 Chinese individuals and replicated them in additional 1,858 Chinese individuals from two cohorts. Over 9% of mQTLs displayed specificity to EAs, facilitating the fine-mapping of EA-specific genetic associations, as shown for variants associated with height. Trans-mQTL hotspots revealed biological pathways contributing to EA-specific genetic associations, including an ERG-mediated 233 trans-mCpG network, implicated in hematopoietic cell differentiation, which likely reflects binding efficiency modulation of the ERG protein complex. More than 90% of mQTLs were shared between different blood cell lineages, with a smaller fraction of lineage-specific mQTLs displaying preferential hypomethylation in the respective lineages. Our study provides new insights into the mQTL landscape across genetic ancestries and their downstream effects on cellular processes and diseases/traits.

A Comprehensive Benchmark of Transcriptomic Biomarkers for Immune Checkpoint Blockades.

Immune checkpoint blockades (ICBs) have revolutionized cancer therapy by inducing durable clinical responses, but only a small percentage of patients can benefit from ICB treatments. Many studies have established various biomarkers to predict ICB responses. However, different biomarkers were found with diverse performances in practice, and a timely and unbiased assessment has yet to be conducted due to the complexity of ICB-related studies and trials. In this study, we manually curated 29 published datasets with matched transcriptome and clinical data from more than 1400 patients, and uniformly preprocessed these datasets for further analyses. In addition, we collected 39 sets of transcriptomic biomarkers, and based on the nature of the corresponding computational methods, we categorized them into the gene-set-like group (with the self-contained design and the competitive design, respectively) and the deconvolution-like group. Next, we investigated the correlations and patterns of these biomarkers and utilized a standardized workflow to systematically evaluate their performance in predicting ICB responses and survival statuses across different datasets, cancer types, antibodies, biopsy times, and combinatory treatments. In our benchmark, most biomarkers showed poor performance in terms of stability and robustness across different datasets. Two scores (TIDE and CYT) had a competitive performance for ICB response prediction, and two others (PASS-ON and EIGS_ssGSEA) showed the best association with clinical outcome. Finally, we developed ICB-Portal to host the datasets, biomarkers, and benchmark results and to implement the computational methods for researchers to test their custom biomarkers. Our work provided valuable resources and a one-stop solution to facilitate ICB-related research.

Dynamic cellular changes in acute kidney injury caused by different ischemia time.

Ischemia reperfusion injury (IRI), often related to surgical procedures, is one of the important causes of acute kidney injury (AKI). To decipher the dynamic process of AKI caused by IRI (with prolonged ischemia phase), we performed single-cell RNA sequencing (scRNA-seq) of clinically relevant IRI murine model with different ischemic intervals. We discovered that Slc5a2hi proximal tubular cells were susceptible to AKI and highly expressed neutral amino acid transporter gene Slc6a19, which was dramatically decreased over the time course. With the usage of mass spectrometry-based metabolomic analysis, we detected that the level of neutral amino acid isoleucine dropped off in AKI mouse plasma metabolites. And the reduction of plasma isoleucine was also verified in patients with cardiac surgery-associated acute kidney injury (CSA-AKI). The findings advanced the understanding of dynamic process of AKI and introduced reduction of isoleucine as a potential biomarker for CSA-AKI.

Rewired m6A epitranscriptomic networks link mutant p53 to neoplastic transformation.

N6-methyladenosine (m6A), one of the most prevalent mRNA modifications in eukaryotes, plays a critical role in modulating both biological and pathological processes. However, it is unknown whether mutant p53 neomorphic oncogenic functions exploit dysregulation of m6A epitranscriptomic networks. Here, we investigate Li-Fraumeni syndrome (LFS)-associated neoplastic transformation driven by mutant p53 in iPSC-derived astrocytes, the cell-of-origin of gliomas. We find that mutant p53 but not wild-type (WT) p53 physically interacts with SVIL to recruit the H3K4me3 methyltransferase MLL1 to activate the expression of m6A reader YTHDF2, culminating in an oncogenic phenotype. Aberrant YTHDF2 upregulation markedly hampers expression of multiple m6A-marked tumor-suppressing transcripts, including CDKN2B and SPOCK2, and induces oncogenic reprogramming. Mutant p53 neoplastic behaviors are significantly impaired by genetic depletion of YTHDF2 or by pharmacological inhibition using MLL1 complex inhibitors. Our study reveals how mutant p53 hijacks epigenetic and epitranscriptomic machinery to initiate gliomagenesis and suggests potential treatment strategies for LFS gliomas.

Benchmark of embedding-based methods for accurate and transferable prediction of drug response.

Prediction of therapy response has been a major challenge in cancer precision medicine due to the extensive tumor heterogeneity. Recently, several deep learning methods have been developed to predict drug response by utilizing various omics data. Most of them train models by using the drug-response screening data generated from cell lines and then use these models to predict response in cancer patient data. In this study, we focus on and evaluate deep learning methods using transcriptome data for the long-standing question of personalized drug-response prediction. We developed an embedding-based approach for drug-response prediction and benchmarked similar methods for their performance. For all methods, we used pretreatment transcriptome data to train models and then conducted a comprehensive evaluation and comparison of the models using cross-panels, cross-datasets and target genes. We further validated the methods using three independent datasets assessing multiple compounds for their predictive capability of drug response, survival outcome and cell line status. As a result, the methods building on gene embeddings had an overall competitive performance with reduced overfitting when we applied evaluation parameters for model fitting as well as the correlation with clinical outcomes in the validation data. We further developed an ensemble model to combine the results from the three most competitive methods for an overall prediction. Finally, we developed DrVAEN (https://bioinfo.uth.edu/drvaen), a user-friendly and easy-accessible web-server that hosts all these methods for drug-response prediction and model comparison for broad use in cancer research, method evaluation and drug development.

Integrated analysis of racial disparities in genomic architecture identifies a trans-ancestry prognostic subtype in bladder cancer.

The incidence of bladder cancer and patient survival vary greatly among different populations, but the influence of the associated molecular features and evolutionary processes on its clinical treatment and prognostication remains unknown. Here, we analyze the genomic architectures of 505 bladder cancer patients from Asian/Black/White populations. We identify a previously unknown association between AHNAK mutations and activity of the APOBEC-a mutational signature, the activity of which varied substantially across populations. All significantly mutated genes but only half of arm-level somatic copy number alterations (SCNAs) are enriched with clonal events, indicating large-scale SCNAs as rich sources of bladder cancer clonal diversities. The prevalence of TP53 and ATM clonal mutations as well as the associated burden of SCNAs is significantly higher in Whites/Blacks than in Asians. We identify a trans-ancestry prognostic subtype of bladder cancer characterized by enrichment of non-muscle-invasive patients and muscle-invasive patients with good prognosis, increased CREBBP/FGFR3/HRAS/NFE2L2 mutations, decreased intra-tumor heterogeneity and genome instability, and an activated tumor microenvironment.

Charting the proteome landscape in major psychiatric disorders: From biomarkers to biological pathways towards drug discovery.

Schizophrenia (SZ), bipolar disorder (BD), and major depressive disorder (MDD) are major mental disorders that affect a significant proportion of the global population. Advancing our knowledge of the pathophysiology of these disorders and identifying biomarkers are urgent needs for developing objective diagnostic tests and new therapeutics. In this study, we performed a systematic review and then extracted, curated, and analyzed proteomics data from published studies, aiming to assess the proteome in peripheral blood of individuals with SZ, BD, or MDD. Then, we performed pathway and network analyses to illuminate the biological themes concatenated by the differentially expressed proteins by systematically interrogating the literature to uncover biological pathways with more robust biological meaning. We identified 486 differentially expressed proteins from 51 studies across the three disorders with 9,423 participants. The great majority of pathways were common to SZ, BD, and MDD. They were related to the immune system, including signaling by interleukins, Toll-like receptor signaling pathway, and complement cascade, and to signal transduction, notably MAPK1/MAPK3 signaling, PI3K-Akt Signaling Pathway, Focal Adhesion-PI3K-Akt-mTOR-signaling, rhodopsin-like receptors, GPCR signaling, and the JAK-STAT signaling pathway. Other shared pathways included advanced glycosylation end-product receptor signaling, Regulation of Insulin-like Growth Factor, cholesterol metabolism, and IL-17 signaling pathway. Pathways shared between SZ and BD were integrin cell-surface interactions, GRB2:SOS provides linkage to MAPK signaling for integrins, and syndecan interactions. Shared between BD and MDD were the NRF2 pathway and signaling by EGFR pathways. Our findings advance our understanding of the protein variations and associations with these disorders, which are useful for accelerating biomarker development and drug discovery.

Hereditary retinoblastoma iPSC model reveals aberrant spliceosome function driving bone malignancies.

The RB1 gene is frequently mutated in human cancers but its role in tumorigenesis remains incompletely defined. Using an induced pluripotent stem cell (iPSC) model of hereditary retinoblastoma (RB), we report that the spliceosome is an up-regulated target responding to oncogenic stress in RB1-mutant cells. By investigating transcriptomes and genome occupancies in RB iPSC­derived osteoblasts (OBs), we discover that both E2F3a, which mediates spliceosomal gene expression, and pRB, which antagonizes E2F3a, coregulate more than one-third of spliceosomal genes by cobinding to their promoters or enhancers. Pharmacological inhibition of the spliceosome in RB1-mutant cells leads to global intron retention, decreased cell proliferation, and impaired tumorigenesis. Tumor specimen studies and genome-wide TCGA (The Cancer Genome Atlas) expression profile analyses support the clinical relevance of pRB and E2F3a in modulating spliceosomal gene expression in multiple cancer types including osteosarcoma (OS). High levels of pRB/E2F3a­regulated spliceosomal genes are associated with poor OS patient survival. Collectively, these findings reveal an undiscovered connection between pRB, E2F3a, the spliceosome, and tumorigenesis, pointing to the spliceosomal machinery as a potentially widespread therapeutic vulnerability of pRB-deficient cancers.

Protein tyrosine phosphatase receptor δ serves as the orexigenic asprosin receptor.

Asprosin is a fasting-induced glucogenic and centrally acting orexigenic hormone. The olfactory receptor Olfr734 is known to be the hepatic receptor for asprosin that mediates its effects on glucose production, but the receptor for asprosin's orexigenic function has been unclear. Here, we have identified protein tyrosine phosphatase receptor δ (Ptprd) as the orexigenic receptor for asprosin. Asprosin functions as a high-affinity Ptprd ligand in hypothalamic AgRP neurons, regulating the activity of this circuit in a cell-autonomous manner. Genetic ablation of Ptprd results in a strong loss of appetite, leanness, and an inability to respond to the orexigenic effects of asprosin. Ablation of Ptprd specifically in AgRP neurons causes resistance to diet-induced obesity. Introduction of the soluble Ptprd ligand-binding domain in the circulation of mice suppresses appetite and blood glucose levels by sequestering plasma asprosin. Identification of Ptprd as the orexigenic asprosin receptor creates a new avenue for the development of anti-obesity therapeutics.

CeDR Atlas: a knowledgebase of cellular drug response.

Drug response to many diseases varies dramatically due to the complex genomics and functional features and contexts. Cellular diversity of human tissues, especially tumors, is one of the major contributing factors to the different drug response in different samples. With the accumulation of single-cell RNA sequencing (scRNA-seq) data, it is now possible to study the drug response to different treatments at the single cell resolution. Here, we present CeDR Atlas (available at https://ngdc.cncb.ac.cn/cedr), a knowledgebase reporting computational inference of cellular drug response for hundreds of cell types from various tissues. We took advantage of the high-throughput profiling of drug-induced gene expression available through the Connectivity Map resource (CMap) as well as hundreds of scRNA-seq data covering cells from a wide variety of organs/tissues, diseases, and conditions. Currently, CeDR maintains the results for more than 582 single cell data objects for human, mouse and cell lines, including about 140 phenotypes and 1250 tissue-cell combination types. All the results can be explored and searched by keywords for drugs, cell types, tissues, diseases, and signature genes. Overall, CeDR fine maps drug response at cellular resolution and sheds lights on the design of combinatorial treatments, drug resistance and even drug side effects.

Cell-Type-Specific Profibrotic Scores across Multi-Organ Systems Predict Cancer Prognosis.

Fibrosis is a major cause of mortality. Key profibrotic mechanisms are common pathways involved in tumorigenesis. Characterizing the profibrotic phenotype will help reveal the underlying mechanisms of early development and progression of a variety of human diseases, such as fibrosis and cancer. Fibroblasts have been center stage in response to various stimuli, such as viral infections. However, a comprehensive catalog of cell types involved in this process is currently lacking. Here, we deployed single-cell transcriptomic data across multi-organ systems (i.e., heart, kidney, liver, and lung) to identify novel profibrotic cell populations based on ECM pathway activity at single-cell resolution. In addition to fibroblasts, we also reported that epithelial, endothelial, myeloid, natural killer T, and secretory cells, as well as proximal convoluted tubule cells of the nephron, were significantly actively involved. Cell-type-specific gene signatures were enriched in viral infection pathways, enhanced glycolysis, and carcinogenesis, among others; they were validated using independent datasets in this study. By projecting the signatures into bulk TCGA tumor samples, we could predict prognosis in the patients using profibrotic scores. Our profibrotic cellular phenotype is useful for identifying new mechanisms and potential drug targets at the cell-type level for a wide range of diseases involved in ECM pathway activation.

An Integrative Transcriptomic and Methylation Approach for Identifying Differentially Expressed Circular RNAs Associated with DNA Methylation Change.

Recently, accumulating evidence has supported that circular RNA (circRNA) plays important roles in tumorigenesis by regulating gene expression at transcriptional and post-transcriptional levels. Expression of circRNAs can be epigenetically silenced by DNA methylation; however, the underlying regulatory mechanisms of circRNAs by DNA methylation remains largely unknown. We explored this regulation in hepatocellular carcinoma (HCC) using genome-wide DNA methylation and RNA sequencing data of the primary tumor and matched adjacent normal tissues from 20 HCC patients. Our pipeline identified 1012 upregulated and 747 downregulated circRNAs (collectively referred to as differentially expressed circRNAs, or DE circRNAs) from HCC RNA-seq data. Among them, 329 DE circRNAs covered differentially methylated sites (adjusted p-value < 0.05, |ΔM| > 0.5) in circRNAs' interior and/or flanking regions. Interestingly, the corresponding parental genes of 46 upregulated and 31 downregulated circRNAs did not show significant expression change in the HCC tumor versus normal samples. Importantly, 34 of the 77 DE circRNAs (44.2%) had significant correlation with DNA methylation change in HCC (Spearman's rank-order correlation, p-value < 0.05), suggesting that aberrant DNA methylation might regulate circular RNA expression in HCC. Our study revealed genome-wide differential circRNA expression in HCC. The significant correlation with DNA methylation change suggested that epigenetic regulation might act on both mRNA and circRNA expression. The specific regulation in HCC and general view in other cancer or disease requires further investigation.

Cell-type deconvolution analysis identifies cancer-associated myofibroblast component as a poor prognostic factor in multiple cancer types.

Cancer-associated fibroblasts (CAFs) constitute a prominent component of the tumor microenvironment and play critical roles in cancer progression and drug resistance. Although recent studies indicate CAFs may consist of several CAF subtypes, the breadth of CAF heterogeneity and functional roles of CAF subtypes in cancer progression remain unclear. In this study, we implemented a cell-type deconvolutional approach to comprehensively characterize cell-type alternations across 18 cancer types from The Cancer Genome Atlas (TCGA). Pan-cancer survival analysis using deconvoluted CAF subtypes revealed myofibroblastic CAF (myCAF) composition as a poor prognostic factor in nine cancer types. Patients with higher myCAF compositions tend to have worse response to six antineoplastic drugs predicted by a lncRNA-based Elastic Net prediction model (LENP). In addition, integrative mutational analysis identified 14 and 413 genes associated with the differentiation degree of myCAF and inflammatory CAF (iCAF), respectively, with significant enrichment of genes involved in fibroblast and extracellular matrix (ECM)-related pathways. In summary, our findings systematically illustrated the complex roles of CAF subtypes in patient prognosis and drug response, and identified putative driver genes in CAF-subtype differentiation. These results provided novel therapeutic perspectives for targeting CAF subtypes in tumor microenvironment and arranging treatment scheme based on the CAF compositions in different cancer types.

DeepVISP: Deep Learning for Virus Site Integration Prediction and Motif Discovery.

Approximately 15% of human cancers are estimated to be attributed to viruses. Virus sequences can be integrated into the host genome, leading to genomic instability and carcinogenesis. Here, a new deep convolutional neural network (CNN) model is developed with attention architecture, namely DeepVISP, for accurately predicting oncogenic virus integration sites (VISs) in the human genome. Using the curated benchmark integration data of three viruses, hepatitis B virus (HBV), human herpesvirus (HPV), and Epstein-Barr virus (EBV), DeepVISP achieves high accuracy and robust performance for all three viruses through automatically learning informative features and essential genomic positions only from the DNA sequences. In comparison, DeepVISP outperforms conventional machine learning methods by 8.43-34.33% measured by area under curve (AUC) value enhancement in three viruses. Moreover, DeepVISP can decode cis-regulatory factors that are potentially involved in virus integration and tumorigenesis, such as HOXB7, IKZF1, and LHX6. These findings are supported by multiple lines of evidence in literature. The clustering analysis of the informative motifs reveales that the representative k-mers in clusters could help guide virus recognition of the host genes. A user-friendly web server is developed for predicting putative oncogenic VISs in the human genome using DeepVISP.

Deep generative neural network for accurate drug response imputation.

Drug response differs substantially in cancer patients due to inter- and intra-tumor heterogeneity. Particularly, transcriptome context, especially tumor microenvironment, has been shown playing a significant role in shaping the actual treatment outcome. In this study, we develop a deep variational autoencoder (VAE) model to compress thousands of genes into latent vectors in a low-dimensional space. We then demonstrate that these encoded vectors could accurately impute drug response, outperform standard signature-gene based approaches, and appropriately control the overfitting problem. We apply rigorous quality assessment and validation, including assessing the impact of cell line lineage, cross-validation, cross-panel evaluation, and application in independent clinical data sets, to warrant the accuracy of the imputed drug response in both cell lines and cancer samples. Specifically, the expression-regulated component (EReX) of the observed drug response achieves high correlation across panels. Using the well-trained models, we impute drug response of The Cancer Genome Atlas data and investigate the features and signatures associated with the imputed drug response, including cell line origins, somatic mutations and tumor mutation burdens, tumor microenvironment, and confounding factors. In summary, our deep learning method and the results are useful for the study of signatures and markers of drug response.

Decoding regulatory structures and features from epigenomics profiles: A Roadmap-ENCODE Variational Auto-Encoder (RE-VAE) model.

The development of chromatin immunoprecipitation (ChIP) with massively parallel DNA sequencing (ChIP-seq) technologies has promoted generation of large-scale epigenomics data, providing us unprecedented opportunities to explore the landscape of epigenomic profiles at scales across both histone marks and tissue types. In addition to many tools directly for data analysis, advanced computational approaches, such as deep learning, have recently become promising to deeply mine the data structures and identify important regulators from complex functional genomics data. We implemented a neural network framework, a Variational Auto-Encoder (VAE) model, to explore the epigenomic data from the Roadmap Epigenomics Project and the Encyclopedia of DNA Elements (ENCODE) project. Our model is applied to 935 reference samples, covering 28 tissues and 12 histone marks. We used the enhancer and promoter regions as the annotation features and ChIP-seq signal values in these regions as the feature values. Through a parameter sweep process, we identified the suitable hyperparameter values and built a VAE model to represent the epigenomics data and to further explore the biological regulation. The resultant Roadmap-ENCODE VAE (RE-VAE) model contained data compression and feature representation. Using the compressed data in the latent space, we found that the majority of histone marks were well clustered but not for tissues or cell types. Tissue or cell specificity was observed only in some histone marks (e.g., H3K4me3 and H3K27ac) and could be characterized when the number of tissue samples is large (e.g., blood and brain). In blood, the contributive regions and genes identified by RE-VAE model were confirmed by tissue-specificity enrichment analysis with an independent tissue expression panel. Finally, we demonstrated that RE-VAE model could detect cancer cell lines with similar epigenomics profiles. In conclusion, we introduced and implemented a VAE model to represent large-scale epigenomics data. The model could be used to explore classifications of histone modifications and tissue/cell specificity and to classify new data with unknown sources.

Molecular signatures identified by integrating gene expression and methylation in non-seminoma and seminoma of testicular germ cell tumours.

Testicular germ cell tumours (TGCTs) are the most common cancer in young male adults (aged 15 to 40). Unlike most other cancer types, identification of molecular signatures in TGCT has rarely reported. In this study, we developed a novel integrative analysis framework to identify co-methylated and co-expressed genes [mRNAs and microRNAs (miRNAs)] modules in two TGCT subtypes: non-seminoma (NSE) and seminoma (SE). We first integrated DNA methylation and mRNA/miRNA expression data and then used a statistical method, CoMEx (Combined score of DNA Methylation and Expression), to assess differentially expressed and methylated (DEM) genes/miRNAs. Next, we identified co-methylation and co-expression modules by applying WGCNA (Weighted Gene Correlation Network Analysis) tool to these DEM genes/miRNAs. The module with the highest average Pearson's Correlation Coefficient (PCC) after considering all pair-wise molecules (genes/miRNAs) included 91 molecules. By integrating both transcription factor and miRNA regulations, we constructed subtype-specific regulatory networks for NSE and SE. We identified four hub miRNAs (miR-182-5p, miR-520b, miR-520c-3p, and miR-7-5p), two hub TFs (MYC and SP1), and two genes (RECK and TERT) in the NSE-specific regulatory network, and two hub miRNAs (miR-182-5p and miR-338-3p), five hub TFs (ETS1, HIF1A, HNF1A, MYC, and SP1), and three hub genes (CDH1, CXCR4, and SNAI1) in the SE-specific regulatory network. miRNA (miR-182-5p) and two TFs (MYC and SP1) were common hubs of NSE and SE. We further examined pathways enriched in these subtype-specific networks. Our study provides a comprehensive view of the molecular signatures and co-regulation in two TGCT subtypes.

KinaseMD: kinase mutations and drug response database.

Mutations in kinases are abundant and critical to study signaling pathways and regulatory roles in human disease, especially in cancer. Somatic mutations in kinase genes can affect drug treatment, both sensitivity and resistance, to clinically used kinase inhibitors. Here, we present a newly constructed database, KinaseMD (kinase mutations and drug response), to structurally and functionally annotate kinase mutations. KinaseMD integrates 679 374 somatic mutations, 251 522 network-rewiring events, and 390 460 drug response records curated from various sources for 547 kinases. We uniquely annotate the mutations and kinase inhibitor response in four types of protein substructures (gatekeeper, A-loop, G-loop and αC-helix) that are linked to kinase inhibitor resistance in literature. In addition, we annotate functional mutations that may rewire kinase regulatory network and report four phosphorylation signals (gain, loss, up-regulation and down-regulation). Overall, KinaseMD provides the most updated information on mutations, unique annotations of drug response especially drug resistance and functional sites of kinases. KinaseMD is accessible at https://bioinfo.uth.edu/kmd/, having functions for searching, browsing and downloading data. To our knowledge, there has been no systematic annotation of these structural mutations linking to kinase inhibitor response. In summary, KinaseMD is a centralized database for kinase mutations and drug response.

A developmental stage-specific network approach for studying dynamic co-regulation of transcription factors and microRNAs during craniofacial development.

Craniofacial development is regulated through dynamic and complex mechanisms that involve various signaling cascades and gene regulations. Disruption of such regulations can result in craniofacial birth defects. Here, we propose the first developmental stage-specific network approach by integrating two crucial regulators, transcription factors (TFs) and microRNAs (miRNAs), to study their co-regulation during craniofacial development. Specifically, we used TFs, miRNAs and non-TF genes to form feed-forward loops (FFLs) using genomic data covering mouse embryonic days E10.5 to E14.5. We identified key novel regulators (TFs Foxm1, Hif1a, Zbtb16, Myog, Myod1 and Tcf7, and miRNAs miR-340-5p and miR-129-5p) and target genes (Col1a1, Sgms2 and Slc8a3) expression of which changed in a developmental stage-dependent manner. We found that the Wnt-FoxO-Hippo pathway (from E10.5 to E11.5), tissue remodeling (from E12.5 to E13.5) and miR-129-5p-mediated Col1a1 regulation (from E10.5 to E14.5) might play crucial roles in craniofacial development. Enrichment analyses further suggested their functions. Our experiments validated the regulatory roles of miR-340-5p and Foxm1 in the Wnt-FoxO-Hippo subnetwork, as well as the role of miR-129-5p in the miR-129-5p-Col1a1 subnetwork. Thus, our study helps understand the comprehensive regulatory mechanisms for craniofacial development.